```{r}
#| message: false
library(readr) # für den data import
library(dplyr) # für das datawrangling
#data_star <-
# readr::read_csv(
# paste0(
# "https://raw.githubusercontent.com/sammerk/",
# "Zuerich_Statistics/master/data/star.csv"
# )
# )
data_star <-
readr::read_csv("data/star.csv")
```Statistische Modellierung mit R
Ein Workshop für die PH Zürich
Methodology: The Science before Statistics
Jede statistische Modellierung gewinnt an Aussagekraft, je umfassender sie die inhaltliche Fragestellung abzubilden im Stande ist. Um aus der riesigen Fülle an Optionen geeignet und zielgerichtet auswählen zu können sind die folgenden Unterscheidungen oft sehr hilfreich.
Erkenntnisinteressen
Ganz grundlegend kann a priori das Erkenntnisinteresse von Studien in die folgenden vier Kategorien unterschieden werden:
| Deskriptiv | Explorativ | Explanativ | Prädiktiv |
|---|---|---|---|
| populationsbeschreibend | hypothesengenerierend | hypothesenprüfend | Datenpunkte vorhersagend oder imputierend |
| Bei welchem Anteil 15-Jähriger in Deutschland handelt es sich um funktionale Analphabet:innen? | Was sind potentielle Ursachen für genderbezogene Disparitäten im Analphabetismus? | Sind 15-jährige Jungen häufiger Analphabeten als 15-jährige Mädchen? | Mit welchen Variablen können Schüler:innen at risk erfolgreich identifiziert werden? |
Gütekriterien wiss. Erkenntnis nach Campbell (1957)
Für ein erfolgreiches Studiendesign und die anschließende statistische Analyse ist es sehr wertvoll sich vorab über Schwerpunkte besonders gewünschter Aspekte wissenschaftlicher Güte Gedanken zu machen. Insbesondere über die Unterkriterien Methodischer Strenge:
- Konstruktvalidität (Inwiefern ist die Interpretation der Messwerte angemessen?)
- Interne Validität (Inwiefern sind Assoziationen von unabhängiger [beeinflussender] und abhängiger [beeinflusster] Variabler als kausale Effekte interpretierbar?)
- Externe Validität (Inwiefern können die Schlussfolgerungen der Studie verallgemeinert werden?)
- Statistische Validität (Wie robust und angemessen sind die verwendeten statistischen Verfahren?)
Deskriptiv- und Inferenzstatistik
Deskriptive Statistik beschreibt vorliegende Daten (z.B. mit Effektstärken), während Inferenzstatistik Aussagen über den die Daten generierenden Mechanismus trifft. Beide können »eher einfach« oder »hoch komplex« sein und oftmals stehen sie in einem synergetischen Verhältnis (siehe Abbildung 1).
Bayesianisches und Frequentistisches Schätzen und Testen
Eine sehr heuristische Klassifikation inferenzstatistischer Verfahren stellt die Unterscheidung von statistischer Schätzung und Testung dar:
(Inferenzstatistische) Schätzungen (estimation with quantified uncertainty) treffen anhand von Stichproben Aussagen über Parameter der Grundgesamtheit (Population) aus der die Stichprobe gezogen wurde.1
1 Bspw.: »Mit 96%iger Wahrscheinlichkeit liegt die Analphabetismusinzidienz von 15-Jährigen in Deutschland zwischen .08 und .12«
(Inferenzstatistische) Hypothesentests bewerten anhand von Stichprobendaten die Gültigkeit von Hypothesen in der Grundgesamtheit (Population) aus der die Stichprobe gezogen wurde.2
2 Bspw.: »Nimmt man an, dass sich die Analphabetismusinzidienz von 15-Jährigen in Deutschland zwischen 2021 und 2025 nicht geändert hat, beträgt die Wahrscheinlichkeit der vorliegenden Daten p = .032«
Diese beiden Verfahren können sowohl im Rahmen der frequentistischen Statistik als auch der bayesianischen Statistik angewendet werden. Die folgende Tabelle gibt einen Überblick über die wichtigsten Werkzeuge:
| Frequentistische Statistik | Bayesianische Statistik | |
|---|---|---|
| Parameterschätzung | Konfidenzintervalle | Posterior Distributions |
| Hypothesentest | p-Werte | Bayes Faktoren & ROPE-CrI Procedure |
Hypothesenarten
Bayesianische wie frequentistischen Hypothesentests können unterschiedliche Arten von Hypothesen zugrunde gelegt werden:
- Punkthypothesen setzen Parameter gleich einer reellen Zahl; etwa \(H_0\text{: } \delta = 0\)
- Äquivalenzhypothesen nehmen Parameter in einem reellen Intervall an; etwa \(H_0\text{: } \delta \not\in\ [-.3, .3]\)
- Informative Hypothesen nehmen eine Ordnungsrelation mehrerer Parameter an; etwa \(\mu_{\text{Baseline}} < \mu_{\text{Imaginary Pill}} < \mu_{\text{Blinded Placebo}}\) (Buergler u. a. 2023)
Die Art der (falsifizierten) Hypothese entscheidet wesentlich stärker über den Informationsgehalt eines Hypothesentests als die Entscheidung für das frequentistische oder bayesianische Paradigma (Hoijtink 2012).
Dies ist am leichtesten anhand der Nullhypothese nachvollziehbar. Wird etwa die Nullhypothese \(H_0\text{: } \delta = 0\) verworfen, wird entsprechend die Alternativhypothese \(H_A\text{: } \delta \neq 0\) angenommen. Diese enthält aber quasi keine Information, da sie nur mit einer einzigen Beobachtung (d = 0.000000 …) verworfen werden kann und im kritischen Rationalismus gilt, dass eine Aussage umso mehr Information enthält, umso leichter sie verworfen werden kann (Döring und Bortz 2016).
Äquivalenzhypothesen können sowohl frequentistisch (z.B. TOAST-Prozedur in R und JASP, Lakens 2017) wie bayesianisch (z.B. ROPE-Ansatz Kruschke 2015) getestet werden. Für das Testen informativer Hypothesen liegen bayesianische Methoden in (u.a.) JASP und R vor (z.B. {bain}, Gu u. a. 2019) sowie in den frequentistischen R-Paketen {restriktor} (Vanbrabant 2020) und {ic.infer} (Grömping 2010).
Grundlagen der Regressionsanalyse
Regressionsanalysen sind ein sehr mächtiges Werkzeug um Zusammenhänge oder Unterschiede zwischen/in Variablen zu modellieren. Innerhalb der Regressionsanalyse kann sowohl geschätzt als auch getestet werden und dies jeweils bayesianisch oder frequentistisch.
Die Grundidee der lineare Regression ist in diesem interaktiven Applet veranschaulicht. Eine Abhängige Variable \(y_i\) wird also in der Regressionsanlyse als normalverteilt mit bedingtem Erwartungswert \(b_0 + b_1*x_i\) und \(\sigma\) angenommen:
\[y_i \sim \mathcal{N}\left( b_0 + b_1*x_i, \sigma \right)\]
Datengrundlage
Das wollen wir in den berühmten Daten aus dem STAR-Experiment (Student/Teacher Achievement Ratio) illustrieren. In diesem Experiment wurden Schüler:innen und Lehrer:innen zufällig auf Klassen mit kleiner (13-17 Schüler:innen pro Lehrer:in) und großer (22-25 Schüler:innen pro Lehrer:in) Klassengröße zugeteilt. Zusätzlich gab es Klassen großer Größe, die von einer:m ausgebildeten Hilfslehrer:in unterstützt wurde. Die Schüler:innen wurden über mehrere Jahre getestet.
Wir können die Daten herunterladen oder direkt programmatisch importieren.
Der Satzsatz sieht wie folgt aus:
```{r}
glimpse(data_star)
```Rows: 26,796
Columns: 18
$ id <dbl> 100017, 100028, 100045, 100045, 100045, 100064, 100070, 100070…
$ sch <dbl> 28, 52, 41, 41, 41, 64, 40, 40, 22, 22, 22, 17, 17, 3, 3, 3, 3…
$ gr <chr> "K", "K", "1", "2", "3", "K", "K", "1", "K", "1", "2", "K", "1…
$ cltype <chr> "small", "reg", "small", "small", "small", "small", "reg", "re…
$ hdeg <chr> "BS/BA", "MS/MA/MEd", "BS/BA", "MS/MA/MEd", "BS/BA", "BS/BA", …
$ clad <chr> "1", "1", "1", "APPR", "1", "1", "APPR", "1", "PROB", "1", "NO…
$ exp <dbl> 3, 12, 20, 15, 5, 19, 2, 5, 9, 9, 2, 20, 16, 8, 23, 21, 14, 6,…
$ trace <chr> "B", "W", "W", "B", "W", "W", "W", "W", "B", "B", "B", "W", "W…
$ read <dbl> 476, 410, 507, 575, 610, 430, 495, 629, 418, 524, 627, 421, 45…
$ math <dbl> 602, 444, 572, 620, 655, 484, 576, 592, 489, 567, 603, 459, 46…
$ ses <chr> "F", "N", "N", "N", "N", "N", "F", "F", "F", "F", "F", "N", "N…
$ schtype <chr> "inner", "suburb", "suburb", "suburb", "suburb", "rural", "rur…
$ sx <chr> "F", "F", "M", "M", "M", "F", "F", "F", "F", "F", "F", "M", "M…
$ eth <chr> "B", "W", "W", "W", "W", "W", "W", "W", "B", "B", "B", "W", "W…
$ birthq <chr> "1980:1", "1980:2", "1980:2", "1980:2", "1980:2", "1980:1", "1…
$ birthy <dbl> 1980, 1980, 1980, 1980, 1980, 1980, 1979, 1979, 1980, 1980, 19…
$ yrs <dbl> 0, 0, 1, 2, 3, 0, 0, 1, 0, 1, 2, 0, 1, 0, 2, 3, 0, 0, 0, 0, 1,…
$ tch <dbl> 478, 893, 698, 701, 706, 1102, 679, 684, 351, 359, 365, 267, 2…Die Variablen dieses Datensatzes sind die folgenden:
id: a factor - student id numbersch: a factor - school id numbergr: grade - an ordered factor with levels K < 1 < 2 < 3cltype: class type - a factor with levels small, reg and reg+A. The last level indicates a regular class size with a teachers aide.hdeg: highest degree obtained by the teacher - an ordered factor with levels ASSOC < BS/BA < MS/MA/MEd < MA+ < Ed.S < Ed.D/Ph.Dclad: career ladder position of the teacher - a factor with levels NOT APPR PROB PEND 1 2 3exp: a numeric vector - the total number of years of experience of the teachertrace: teacher’s race - a factor with levels W, B, A, H, I and O representing white, black, Asian, Hispanic, Indian (Native American) and otherread: the student’s total reading scaled scoremath: the student’s total math scaled scoreses: socioeconomic status - a factor with levels F and N representing eligible for free lunches or not eligibleschtype: school type - a factor with levels inner, suburb, rural and urbansx: student’s sex - a factor with levels M Feth: student’s ethnicity - a factor with the same levels as tracebirthq: student’s birth quarter - an ordered factor with levels 1977:1 < … < 1982:2birthy: student’s birth year - an ordered factor with levels 1977:1982yrs: number of years of schooling for the student - a numeric version of the grade gr with Kindergarten represented as 0. This variable was generated from gr and does not allow for a student being retained.tch: a factor - teacher id number
Um nicht gleich mit Mehrebenenregression starten zu müssen, erstellen wir als erstes einen Subdatensatz, der nur eine:n Schüler:in pro Lehrer:in enthält
```{r}
data_star_sub <-
data_star %>%
group_by(tch) %>%
sample_n(1) %>%
ungroup()
```»Effekte« der Klassenstufe
Einfache lineare Regression
Zunächst wollen wir »Effekte« der Klassenstufe modellieren. Ich empfehle dingend Daten immer erst zu visualisieren und dann statistisch zu modellieren.
```{r}
library(ggplot2) # plots
```Warning: package 'ggplot2' was built under R version 4.4.1```{r}
library(ggforce) # sina plots
ggplot(data_star_sub, aes(gr, math)) +
geom_jitter() +
theme_minimal()
```Warning: Removed 112 rows containing missing values or values outside the scale range
(`geom_point()`).
Die unabhängige Variable ist hier noch nicht intervallskaliert. Dies können wir wie folgt ändern.
```{r}
# Das Schuljahr K zu 0 rekodieren
data_star_sub <-
data_star_sub %>%
mutate(
grade = case_when(
gr == "K" ~ 0,
gr == "1" ~ 1,
gr == "2" ~ 2,
gr == "3" ~ 3
)
)
# Anschließend nochmals visualisieren
ggplot(data_star_sub, aes(grade, math)) +
geom_jitter() +
theme_minimal()
```Warning: Removed 112 rows containing missing values or values outside the scale range
(`geom_point()`).
Dann kann eine erste Regressionsgerade modelliert werden.
Frequentistische Punktschätzung
Mit der Funktion lm() werden mithilfe kleinster Quadrate \(b_0\) und \(b_1\) geschätzt.
```{r}
mod0 <- lm(math ~ grade, data = data_star_sub)
mod0
```
Call:
lm(formula = math ~ grade, data = data_star_sub)
Coefficients:
(Intercept) grade
485.08 45.15 Frequentistischer Nullhypothesensignifikanztest
Die Standardinferenzstatistik zu einem solchen Modell ist sind p-Werte für die Nullhypothesen \[H_0:\; b_0 = 0\] \[H_0:\; b_1 = 0\]
```{r}
summary(mod0)
```
Call:
lm(formula = math ~ grade, data = data_star_sub)
Residuals:
Min 1Q Median 3Q Max
-114.515 -29.869 -4.515 26.777 145.777
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 485.077 2.056 235.92 <2e-16 ***
grade 45.146 1.126 40.08 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 44.13 on 1273 degrees of freedom
(112 observations deleted due to missingness)
Multiple R-squared: 0.5579, Adjusted R-squared: 0.5575
F-statistic: 1606 on 1 and 1273 DF, p-value: < 2.2e-16Frequentistische Konfidenzintervalle
```{r}
confint(mod0)
``` 2.5 % 97.5 %
(Intercept) 481.04329 489.11078
grade 42.93623 47.35601Bayesianische Posterior Distributions
```{r}
#| message: false
#| warning: false
library(MCMCpack) #für bayesianische Schätzung
mod1 <- MCMCregress(math ~ grade, data = data_star_sub)
plot(mod1)
summary(mod1)
```
Iterations = 1001:11000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 10000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
(Intercept) 485.08 2.056 0.02056 0.02056
grade 45.14 1.139 0.01139 0.01139
sigma2 1950.64 77.497 0.77497 0.77497
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
(Intercept) 481.05 483.69 485.10 486.47 489.12
grade 42.92 44.37 45.13 45.89 47.38
sigma2 1804.07 1897.29 1948.67 2002.21 2106.36
Während die Densityplots den Posterior der Schätzung charakterisieren, erhält man mit summary() auch die 95%CrI-Intervalle.
Bayes Factor
Mit einem Bayes Factor testet man im Falle unseres Models wie wahrscheinlich die vorliegenden Daten sind geben ein Intercept-Only-Modell vs. unser Model mod1. Da die Funktion lmBF() im Paket {BayesFactor} allerdings Daten ohne Missings erwartet, müssen diese noch gefiltert werden.
```{r}
#| message: false
library(BayesFactor)
```Warning: package 'Matrix' was built under R version 4.4.1```{r}
#| message: false
mod2 <- lmBF(
formula = math ~ grade,
data = data_star_sub %>%
filter(!is.na(math))
)
```Warning: data coerced from tibble to data frameStandardisierte Regression
Die in mod0 erhaltenen Regressionskoeffizienten konnten wir nutzen um zu interpretieren welche Punktzahl durchschnittlich im Kindergarten vorliegt (\(b_0\)) und was der durchschnittliche Anstieg in einem Jahr ist (\(b_1\)). Diese Zahlen sind insbesondere dann gut als (nicht-standardisierte) Effektstärke interpretierbar, wenn man mit der Skala vertraut ist. Eine Alternative dazu ist die z-Standardisierung der abhängigen Variable. Um dies zu veranschaulichen betrachten wir Effekte des SES (free lunch) auf die Mathematikleistung in der dritten Klasse. Die Rohdataten sehen wie folgt aus:
```{r}
data_star_sub %>%
filter(grade == 3 & !is.na(ses)) %>%
ggplot(aes(ses, math)) +
geom_violin() +
geom_sina() +
stat_summary(
fun.data = mean_sdl,
fun.args = list(mult = 1),
geom = "pointrange",
colour = "#1bbc9d"
) +
theme_minimal()
```Warning: Removed 36 rows containing non-finite outside the scale range
(`stat_ydensity()`).Warning: Removed 36 rows containing non-finite outside the scale range
(`stat_sina()`).Warning: Removed 36 rows containing non-finite outside the scale range
(`stat_summary()`).
z-standardisieren wir die abhängige Variable bleiben die Verteilungsformen gleich, lediglich die Skalierung ändert sich.
```{r}
data_star_sub %>%
filter(grade == 3 & !is.na(ses)) %>%
ggplot(aes(ses, scale(math))) +
geom_violin() +
geom_sina() +
stat_summary(
fun.data = mean_sdl,
fun.args = list(mult = 1),
geom = "pointrange",
colour = "#1bbc9d"
) +
theme_minimal()
```Warning: Removed 36 rows containing non-finite outside the scale range
(`stat_ydensity()`).Warning: Removed 36 rows containing non-finite outside the scale range
(`stat_sina()`).Warning: Removed 36 rows containing non-finite outside the scale range
(`stat_summary()`).
Der folgende Code kodiert automatisch die ses Variable als 0, 1-Variable um.
```{r}
mod3 <-
lm(
scale(math) ~ ses,
data = data_star_sub %>%
filter(grade == 3 & !is.na(ses))
)
```Die Koeffizienten passend sehr gut zur Abbildung:
\[ \operatorname{\widehat{scale(math)}} = -0.4 + 0.76(\operatorname{ses}_{\operatorname{N}}) \]
-0.4 ist der Mittelwert der Gruppe mit Free Lunch, 0.76 die Differenz der beiden Gruppenmittelwerte gemessen in Standardabweichungen und damit wie ein Cohen’s \(d\) interpretiertbar.
```{r}
library(effectsize)
```Warning: package 'effectsize' was built under R version 4.4.1```{r}
cohens_d(
math ~ ses,
data = data_star_sub %>% filter(grade == 3 & !is.na(ses)),
pooled_sd = F
)
```Warning: Missing values detected. NAs dropped.Cohen's d | 95% CI
--------------------------
-0.82 | [-1.06, -0.58]
- Estimated using un-pooled SD.Multiple Regression
Die multiple Regression erweitert die Idee der einfachen lineare Regression auf mehrere unabhängige Variablen:
\[y_i \sim \mathcal{N}\left( b_0 + b_1*x_{1i} + b_2*x_{2i}, \sigma \right)\]
Literatur
Buergler, Sarah, Dilan Sezer, Niels Bagge, Irving Kirsch, Cosima Locher, Claudia Carvalho, und Jens Gaab. 2023. „Imaginary Pills and Open-Label Placebos Can Reduce Test Anxiety by Means of Placebo Mechanisms Scientific Reports“. Scientific Reports 13 (1): 2624. https://doi.org/10.1038/s41598-023-29624-7.
Campbell, Donald T. 1957. „Factors Relevant to the Validity of Experiments in Social Settings“. Psychological Bulletin 54 (4): 297–312. https://doi.org/10.1037/h0040950.
Döring, Nicola, und Jürgen Bortz. 2016. Forschungsmethoden und Evaluation in den Sozial- und Humanwissenschaften. 5., vollst. Berlin, Heidelberg: Springer.
Grömping, Ulrike. 2010. „Inference with Linear Equality and Inequality Constraints Using R: The Package Ic.infer“. Journal of Statistical Software 33 (Februar): 1–31. https://doi.org/10.18637/jss.v033.i10.
Gu, Xin, Herbert Hoijtink, Joris Mulder, und Yves Rosseel. 2019. „Bain: A Program for Bayesian Testing of Order Constrained Hypotheses in Structural Equation Models“. Journal of Statistical Computation and Simulation 89 (8): 1526–53. https://doi.org/10.1080/00949655.2019.1590574.
Hoijtink, Herbert. 2012. Informative Hypotheses: Theory and Practice for Behavioral and Social Scientists. Chapman & Hall/CRC Statistics in the Social and Behavioral Sciences Series. Boca Raton: CRC.
Kruschke, John K. 2015. Doing Bayesian Data Analysis: A Tutorial with R, JAGS, and Stan. 2nd Aufl. Academic Press.
Lakens, Daniël. 2017. „Equivalence Tests: A Practical Primer for t Tests, Correlations, and Meta-Analyses“. Social Psychological and Personality Science 8 (4): 355–62. https://doi.org/10.1177/1948550617697177.
Vanbrabant, L. 2020. Restriktor: Constrained Statistical Inference. https://cran.r-project.org/package=restriktor.